#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#include <time.h>

#define G 2.959572553206070e-004
#define diretorioWindows "C:\\dados\\"

void criarMatrizRK45(void);
void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho);
void abrirArquivos(char SO,char* metodo,double tFinal,double n);
void carregarArquivosLinha(double *z, int dimensao,int tamanho);
void FD(double *J,double* delta,int tamanho,double *resultado);
void FZ(double *dx, double *z, double *distancia, double *m, int dimensao);
void produtoEscalarVetor(double* v,double a,int tamanho, double *resultado);
void somaVetores(double* v,double* v2,int tamanho, double *resultado);
void imprimir(double* v,int tamanho);
double calculaNorma(double *v,int dimensao);
void calcularDistancia(double *z,double *distancia, int dimensao, int corpos);
void finalizar(void);

FILE *f[7];
double *z,*J,*RD1,*RD2,*RD3,*RD4,*RD5,*RD6,*RZ1,*RZ2,*RZ3,*RZ4,*RZ5,*RZ6,*delta,*temp,*temp2,temp4,*norma,*Y,*y,*distancia,*m,*a,*b,*c2;
int corpos;

int main(int argc, char *argv[])
{
	clock_t inicio=clock();

	double tFinal,n,dt;
	int i=0,j,dimensao,T,tamanho,renormalizacao,escrever;
	char aux[100];

	escrever=atoi(argv[1]);
	tFinal=atof(argv[2]);
	n=atof(argv[3]);

	dimensao=2;
	corpos=3;
	tamanho=2*corpos*dimensao;
	T=pow(2,n)+1;
	renormalizacao=1;
	dt=tFinal/pow(2,n);
	criarMatrizRK45();

	abrirArquivos('w',"artigo2",tFinal,n);
	J=(double *)calloc(pow(tamanho,2), sizeof(double));
	delta=(double *)calloc(tamanho, sizeof(double));
	temp=(double *)calloc(tamanho, sizeof(double));
	temp2=(double *)calloc(tamanho, sizeof(double));
	norma=(double*)calloc(2,sizeof(double));
	distancia=(double*)calloc(corpos*corpos,sizeof(double));
	Y=(double*)calloc(2,sizeof(double));
	y=(double*)calloc(2,sizeof(double));
	z=(double*)calloc(tamanho,sizeof(double));
	m=(double*)calloc(corpos,sizeof(double));
	RD1=(double*)calloc(tamanho,sizeof(double));
	RD2=(double*)calloc(tamanho,sizeof(double));
	RD3=(double*)calloc(tamanho,sizeof(double));
	RD4=(double*)calloc(tamanho,sizeof(double));
	RD5=(double*)calloc(tamanho,sizeof(double));
	RD6=(double*)calloc(tamanho,sizeof(double));
	RZ1=(double*)calloc(tamanho,sizeof(double));
	RZ2=(double*)calloc(tamanho,sizeof(double));
	RZ3=(double*)calloc(tamanho,sizeof(double));
	RZ4=(double*)calloc(tamanho,sizeof(double));
	RZ5=(double*)calloc(tamanho,sizeof(double));
	RZ6=(double*)calloc(tamanho,sizeof(double));

	m[0]=0.3;
	m[1]=9.543512141169373e-004;
	m[2]=3.003167261575587e-006;

	sprintf(aux,"%0.15lf ",Y[0]);
	fputs(aux,f[6]);
	fputs("\n",f[6]);

	//delta inicial aleatorio
	srand( (unsigned)time(NULL) );
	for(i=0;i<tamanho;i++)
	{
		//delta[i]=rand()/(double)RAND_MAX;
		//delta[i]=delta[i]-0.5;
		//delta[i]=delta[i]*0.01;
		delta[i]=1.0;
	}
	//delta[4]=0.1;

	// norma delta inicial nao normalizado
	norma[0]=calculaNorma(delta,tamanho);
	//normalizacao delta inicial
	for(i=0;i<tamanho;i++)
	{
		delta[i]=delta[i]/norma[0];
	}
	imprimir(delta,12);

	//printf("\n%d\n",renormalizacao);

	//norma delta inicial normalizado
	norma[0]=calculaNorma(delta,tamanho);

	puts("");

	for(i=1;i<T;i++)
	{
		carregarArquivosLinha(z,dimensao,tamanho);

		//stage 1
		FZ(temp, z, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ1);

		jacobiana3C(J,m[0],m[1],m[2],z[0],z[1],z[2],z[3],z[4],z[5],tamanho);
		FD(J,delta,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD1);

		//stage 2
		produtoEscalarVetor(RZ1,a[1*6+0],tamanho,temp);
		//temp2- z intermediario
		somaVetores(temp,z,tamanho,temp2);
		FZ(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ2);

		jacobiana3C(J,m[0],m[1],m[2],temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],tamanho);
		produtoEscalarVetor(RD1,a[1*6+0],tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		FD(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD2);

		//stage 3
		produtoEscalarVetor(RZ1,a[2*6+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(RZ2,a[2*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);

		FZ(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ3);

		jacobiana3C(J,m[0],m[1],m[2],temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],tamanho);
		produtoEscalarVetor(RD1,a[2*6+0],tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(RD2,a[2*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FD(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD3);

		//stage 4
		produtoEscalarVetor(RZ1,a[3*6+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(RZ2,a[3*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ3,a[3*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FZ(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ4);

		jacobiana3C(J,m[0],m[1],m[2],temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],tamanho);
		produtoEscalarVetor(RD1,a[3*6+0],tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(RD2,a[3*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD3,a[3*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FD(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD4);

		//stage 5
		produtoEscalarVetor(RZ1,a[4*6+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(RZ2,a[4*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ3,a[4*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ4,a[4*6+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FZ(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ5);

		jacobiana3C(J,m[0],m[1],m[2],temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],tamanho);
		produtoEscalarVetor(RD1,a[4*6+0],tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(RD2,a[4*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD3,a[4*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD4,a[4*6+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FD(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD5);

		//stage 6
		produtoEscalarVetor(RZ1,a[5*6+0],tamanho,temp);
		somaVetores(temp,z,tamanho,temp2);
		produtoEscalarVetor(RZ2,a[5*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ3,a[5*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ4,a[5*6+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RZ5,a[5*6+4],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FZ(temp, temp2, distancia, m, dimensao);
		produtoEscalarVetor(temp,dt,tamanho,RZ6);

		jacobiana3C(J,m[0],m[1],m[2],temp2[0],temp2[1],temp2[2],temp2[3],temp2[4],temp2[5],tamanho);
		produtoEscalarVetor(RD1,a[5*6+0],tamanho,temp);
		somaVetores(temp,delta,tamanho,temp2);
		produtoEscalarVetor(RD2,a[5*6+1],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD3,a[5*6+2],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD4,a[5*6+3],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		produtoEscalarVetor(RD5,a[5*6+4],tamanho,temp);
		somaVetores(temp,temp2,tamanho,temp2);
		FD(J,temp2,tamanho,temp);
		produtoEscalarVetor(temp,dt,tamanho,RD6);

		produtoEscalarVetor(RD1,b[0],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);
		produtoEscalarVetor(RD2,b[1],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);
		produtoEscalarVetor(RD3,b[2],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);
		produtoEscalarVetor(RD4,b[3],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);
		produtoEscalarVetor(RD5,b[4],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);
		produtoEscalarVetor(RD6,b[5],tamanho,temp);
		somaVetores(delta,temp,tamanho,delta);

		norma[1]=calculaNorma(delta,tamanho);

		y[1]=((double)(i-1)/(double)(i))*y[0]+2*log(norma[1]/norma[0]);
		Y[1]=((double)(i-1)*Y[0]+y[1])/(double)i;

		if(i%renormalizacao==0)
		{
			for(j=0;j<tamanho;j++)
			{
				delta[j]=delta[j]/norma[1];
			}
			norma[1]=1;
		}

		if(i%escrever==0)
		{
			sprintf(aux,"%0.15lf ",Y[1]);
			fputs(aux,f[6]);
			fputs("\n",f[6]);
		}

		y[0]=y[1];
		Y[0]=Y[1];
		norma[0]=norma[1];

	}
	printf("%0.20lf\n",Y[1]);

	int minutos=0;
	int segundos=0;
	int horas=0;
	double tempo=0;
	int tempo2=0;

	clock_t fim=clock();
	clock_t total=fim-inicio;
	tempo=(double)total/CLOCKS_PER_SEC;
	tempo2=(int)tempo;
	horas=(tempo2/3600);
	minutos=(tempo2/60)-horas*60;
	segundos=(tempo2 % 60);
	printf("%d h %d min %d s\n",horas, minutos ,segundos);
	printf("%lf \n",tempo);

	void finalizar(void);
	return 0;
}

void imprimir(double* v,int tamanho)
{
	int i;
	for(i=0;i<tamanho;i++)
	{
		printf("%0.10lf\n",v[i]);
	}
	puts("");
}

void somaVetores(double* v,double* v2,int tamanho, double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=v[i]+v2[i];
	}
}

void produtoEscalarVetor(double* v,double a,int tamanho,double *resultado)
{
	int i;

	for(i=0;i<tamanho;i++)
	{
		resultado[i]=a*v[i];
	}
}

void FZ(double *dx, double *z, double *distancia, double *m, int dimensao)
{

	dx[0]=z[6]/m[0];
	dx[1]=z[7]/m[0];
	dx[2]=z[8]/m[1];
	dx[3]=z[9]/m[1];
	dx[4]=z[10]/m[2];
	dx[5]=z[11]/m[2];

	calcularDistancia(z,distancia,dimensao,corpos);

	dx[6]=      (G*m[0]*m[1]*(z[2]-z[0]))/(double)(pow(distancia[0*corpos+1],3));
	dx[6]=dx[6]+(G*m[0]*m[2]*(z[4]-z[0]))/(double)(pow(distancia[0*corpos+2],3));
	dx[7]=      (G*m[0]*m[1]*(z[3]-z[1]))/(double)(pow(distancia[0*corpos+1],3));
	dx[7]=dx[7]+(G*m[0]*m[2]*(z[5]-z[1]))/(double)(pow(distancia[0*corpos+2],3));

	dx[8]=      (G*m[1]*m[0]*(z[0]-z[2]))/(double)(pow(distancia[1*corpos+0],3));
	dx[8]=dx[8]+(G*m[1]*m[2]*(z[4]-z[2]))/(double)(pow(distancia[1*corpos+2],3));
	dx[9]=      (G*m[1]*m[0]*(z[1]-z[3]))/(double)(pow(distancia[1*corpos+0],3));
	dx[9]=dx[9]+(G*m[1]*m[2]*(z[5]-z[3]))/(double)(pow(distancia[1*corpos+2],3));

	dx[10]=       (G*m[2]*m[0]*(z[0]-z[4]))/(double)(pow(distancia[2*corpos+0],3));
	dx[10]=dx[10]+(G*m[2]*m[1]*(z[2]-z[4]))/(double)(pow(distancia[2*corpos+1],3));
	dx[11]=       (G*m[2]*m[0]*(z[1]-z[5]))/(double)(pow(distancia[2*corpos+0],3));
	dx[11]=dx[11]+(G*m[2]*m[1]*(z[3]-z[5]))/(double)(pow(distancia[2*corpos+1],3));

}

void FD(double *J,double* delta,int tamanho, double *resultado)
{
	int i,j;
	double temp;

	for(i=0;i<tamanho;i++)
	{
		temp=J[i*tamanho]*delta[0];
		for(j=1;j<tamanho;j++)
		{
			temp+=J[i*tamanho+j]*delta[j];
		}
		resultado[i]=temp;
	}
}

void abrirArquivos(char SO,char* metodo,double tFinal,double n)
{
	char *barra,*diretorio,*temp,*raiz;

	raiz=(char*)calloc(1000,sizeof(char));
	diretorio=(char*)calloc(1000,sizeof(char));
	temp=(char*)calloc(1000,sizeof(char));

	if(SO=='w')
	{
		barra="\\";
	}
	else
	{
		barra="/";
	}

	strcpy(raiz,diretorioWindows);
	strcat(raiz,metodo);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", tFinal);
	strcat(raiz,temp);
	strcat(raiz,barra);
	sprintf(temp, "%0.0lf", n);
	strcat(raiz,temp);
	strcat(raiz,barra);

	strcpy(diretorio,raiz);
	strcat(diretorio,"r1.txt");
	f[0]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"r2.txt");
	f[1]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"r3.txt");
	f[2]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p1.txt");
	f[3]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p2.txt");
	f[4]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"p3.txt");
	f[5]=fopen(diretorio,"r");

	strcpy(diretorio,raiz);
	strcat(diretorio,"megnoRK5.txt");
	f[6]=fopen(diretorio,"w");

	free(raiz);
	free(diretorio);
	free(temp);
}

void carregarArquivosLinha(double *z, int dimensao,int tamanho)
{
	int k=0;
	char *temp,*temp2;

	for(k=0;k<tamanho/2;k++)
	{
		temp=(char*)calloc(1000,sizeof(char));
		if(!temp)
		{
			puts("carregar temp");
		}
		temp2=(char*)calloc(1000,sizeof(char));
		if(!temp2)
		{
			puts("carregar temp2");
		}
		fgets(temp,1000,f[k]);
		int i=0,j=0;
		do{
			temp2[i]=temp[i];
			i++;
		}while(temp[i]!=' ');
		z[2*k]=atof(temp2);
		do{
			temp2[j]=temp[i];
			i++;
			j++;
		}while(temp[i]!=' ');
		z[2*k+1]=atof(temp2);
		free(temp);
		free(temp2);
	}
}

void jacobiana3C(double* j, double m1,double m2,double m3,double r1x,double r1y,double r2x,double r2y,double r3x,double r3y,int tamanho)
{
	double G2=(double)G;
	double t70 = r1x-r2x;
	double t71 = r1y-r2y;
	double t72 = r1x-r3x;
	double t73 = r1y-r3y;
	double t78 = r1x*2;
	double t83 = r2x*2;
	double t74 = t78-t83;
	double t75 = pow(t70,2);
	double t76 = pow(t71,2);
	double t77 = t75+t76;
	double t85 = r3x*2;
	double t79 = -1*t78+t85;
	double t80 = pow(t72,2);
	double t81 = pow(t73,2);
	double t82 = t80+t81;
	double t84 = pow((1/t77),(5.0/2.0));
	double t86 = r1y*2;
	double t87 = pow((1/t82),(5.0/2.0));
	double t88 = pow((1/t77),(3.0/2.0));
	double t89 = pow(t74,2);
	double t90 = G2 * m1 * m2 * t84 * t89* 3.0/4.0;
	double t91 = r2y*2;
	double t92 = t86-t91;
	double t93 = G2*m1*m2*t74*t84*t92*3.0/4.0;
	double t94 = pow((1.0/t82),(3.0/2.0));
	double t95 = t78-t85;
	double t96 = r3y*2.0;
	double t97 = t86-t96;
	double t98 = G2*m1*m3*t87*t95*t97*(3.0/4.0);
	double t99 = t93+t98;
	double t100 = G2*m1*m2*t88;
	double t101 = G2*m1*m3*t94;
	double t102 = pow(t92,2);
	double t103 = G2*m1*m2*t102*t84*(3.0/4.0);
	double t104 = pow(t97,2);
	double t105 = G2*m1*m3*t104*t87*(3.0/4.0);
	double t106 = t100-t90;
	double t107 = r2x-r3x;
	double t108 = r2y-r3y;
	double t109 = t83-t85;
	double t110 = pow(t107,2);
	double t111 = pow(t108,2);
	double t112 = t110+t111;
	double t113 = pow((1.0/t112),(5.0/2.0));
	double t114 = pow((1.0/t112),(3.0/2.0));
	double t115 = pow(t109,2);
	double t116 = G2*m2*m3*t113*t115*(3.0/4.0);
	double t117 = t91-t96;
	double t118 = G2*m2*m3*t109*t113*t117*(3.0/4.0);
	double t119 = t100-t103;
	double t120 = t118+t93;
	double t121 = G2*m2*m3*t114;
	double t122 = pow(t117,2);
	double t123 = G2*m2*m3*t113*t122*(3.0/4.0);
	double t124 = pow(t95,2);
	double t127 = G2*m1*m3*t124*t87*(3.0/4.0);
	double t125 = t101-t127;
	double t126 = -1*t116+t121;
	double t128 = t101-t105;
	double t129 = t121-t123;
	double t130 = t118+t98;

	int colunas=3*2*2;

	j[0*colunas+6]=1/m1;
	j[1*colunas+7]=1/m1;
	j[2*colunas+8]=1/m2;
	j[3*colunas+9]=1/m2;
	j[4*colunas+10]=1/m3;
	j[5*colunas+11]=1/m3;
	j[6*colunas+0]=t90-G2*m1*m2*t88-G2*m1*m3*t94+G2*m1*m3*pow(t79,2)*t87*(3.0/4.0);
	j[6*colunas+1]=t99;
	j[6*colunas+2]=t106;
	j[6*colunas+3]=-t93;
	j[6*colunas+4]=t125;
	j[6*colunas+5]=-t98;
	j[7*colunas+0]=j[0*colunas+1];
	j[7*colunas+1]=-t100-t101+t103+t105;
	j[7*colunas+2]=-t93;
	j[7*colunas+3]=t119;
	j[7*colunas+4]=-t98;
	j[7*colunas+5]=t128;
	j[8*colunas+0]= t106;
	j[8*colunas+1]=-t93;
	j[8*colunas+2]=-t100+t116+t90-G2*m2*m3*t114;
	j[8*colunas+3]=t120;
	j[8*colunas+4]=t126;
	j[8*colunas+5]=-t118;
	j[9*colunas+0]=-t93;
	j[9*colunas+1]=t119;
	j[9*colunas+2]=t120;
	j[9*colunas+3]=-t100+t103-t121+t123;
	j[9*colunas+4]=-t118;
	j[9*colunas+5]=t129;
	j[10*colunas+0]=t125;
	j[10*colunas+1]=-t98;
	j[10*colunas+2]=t126;
	j[10*colunas+3]=-t118;
	j[10*colunas+4]=-t101+t116-t121+t127;
	j[10*colunas+5]=t130;
	j[11*colunas+0]=-t98;
	j[11*colunas+1]=t128;
	j[11*colunas+2]=-t118;
	j[11*colunas+3]=t129;
	j[11*colunas+4]=t130;
	j[11*colunas+5]=-t101+t105-t121+t123;

}

double calculaNorma(double *v,int dimensao)
{
	double temp=0;
	int i;

	for(i=0;i<dimensao;i++)
	{
		temp+=v[i]*v[i];
	}
	temp=sqrt(temp);
	return temp;
}

void calcularDistancia(double *z,double *distancia, int dimensao, int corpos)
{
	int i,j,k;

	for(i=0;i<corpos;i++)
	{
		for(j=0;j<corpos;j++)
		{
			distancia[i*corpos+j]=0;
			for(k=0;k<dimensao;k++)
			{
				temp4=(z[i*dimensao+k]-z[j*dimensao+k]);
				temp4=temp4*temp4;
				distancia[i*corpos+j]=distancia[i*corpos+j]+temp4;
			}
			distancia[i*corpos+j]=(double)sqrt((double)distancia[i*corpos+j]);
		}
	}
}

void criarMatrizRK45(void)
{
	int stages=6;
	a=(double*)calloc(stages*stages, sizeof(double));
	b=(double*)calloc(stages, sizeof(double));
	c2=(double*)calloc(stages, sizeof(double));

	a[1*stages+0]=1.0/4;
	a[2*stages+0]=3.0/32;
	a[2*stages+1]=9.0/32;
	a[3*stages+0]=1932.0/2197;
	a[3*stages+1]=-7200.0/2197;
	a[3*stages+2]=7296.0/2197;
	a[4*stages+0]=439.0/216;
	a[4*stages+1]=-8.0;
	a[4*stages+2]=3680.0/513;
	a[4*stages+3]=-845.0/4104;
	a[5*stages+0]=-8.0/27;
	a[5*stages+1]=2.0;
	a[5*stages+2]=-3544.0/2565;
	a[5*stages+3]=1859.0/4104;
	a[5*stages+4]=-11.0/40;

	//final
	b[0]=16.0/135;
	b[1]=0.0;
	b[2]=6656.0/12825;
	b[3]=28561.0/56430;
	b[4]=-9.0/50;
	b[5]=2.0/55;

	//independente
	c2[0]=0.0;
	c2[1]=1.0/4;
	c2[2]=3.0/8;
	c2[3]=12.0/13;
	c2[4]=1.0;
	c2[5]=1.0/2;
}

void finalizar(void)
{
	fclose(f[0]);
	fclose(f[1]);
	fclose(f[2]);
	fclose(f[3]);
	fclose(f[4]);
	fclose(f[5]);
	fclose(f[6]);

	free(distancia);
	free(m);
	free(a);
	free(b);
	free(c2);
	free(RZ1);
	free(RZ2);
	free(RZ3);
	free(RZ4);
	free(RZ5);
	free(RZ6);
	free(RD1);
	free(RD2);
	free(RD3);
	free(RD4);
	free(RD5);
	free(RD6);
	free(J);
	free(delta);
	free(temp);
	free(temp2);
	free(norma);
	free(Y);
	free(y);
	free(z);
}


